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Abstract 

The  problem  of  estimating  Lhe  orbit  parameters  from  oarly-orbit 
observations  of  an  earth  satellite  is  used  to  compare  the  accuracy  and 
application  of  the  Extended  Kalman  filter  and  the  classical  filtering 
method  of  Weighted  Least  Squares.  To  obtain  an  absolute  comparison, 
a  true  two-body,  drag-free  Keplerian  orbit  is  simulated,  observations 
are  computed  and  contaminated  with  noise,  and  the  orbit  parameters 
estimated  by  each  filter  are  compared.  The  accuracy  of  the  two  filters 
was  compared  using  the  same  set  of  observations  to  determine  the  effect 
of  observation  truncation  and  initial  conditions  on  the  results. 

Based  on  this  study  it  was  concluded  that  the  Weighted  Least  Squares 
filter  and  the  Extended  Kalman  filter  yield  about  the  same  accuracy  in 
the  oarly-orbit  determination  problem. 
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A  COMPARISON  OP  THE  EXTENDED  KALMAN 
FILTER  AND  WEIGHTED  LEAST  SQUARES 
IN  EARLY-ORBIT  DETERMINATION 

I .  Introduction 

The  determination  of  the  orbit  parameters  of  man-made  satellites 
usually  involves  the  calculation  of  the  satellite's  position  and 
velocity  at  some  epoch  time  using  a  number  of  observations  from  the 
earth's  surface.  These  observations  are  inherently  noisy,  so  the 
determination  of  the  orbit  parameters  is  a  statistical  problem,  thus 
a  "best"  estimate  of  the  orbit  parameters  is  obtained  using  some  form 
of  a  statistical  filter. 

One  of  the  traditional  methods  of  solving  this  filter  problem  is 
via  differential  correction  of  initial  conditions  using  the  classical 
method  of  Weighted  Least  Squares  (WLS).  This  method  produces  a  very 
good  estimate  of  the  orbit  parameters  when  a  large  number  of  observa¬ 
tions  are  available,  i.c.,  when  the  system  is  highly  overdetermined. 
Problems  arise  when  the  number  of  observations  is  small,  thus  one  of 
the  goals  of  this  study  is  to  investigate  the  application  of  the  WLS 
filter  to  early-orbit  determination.  Early- orbit  determination  is 
extremely  important  as  the  orbit  parameters  must  be  quickly  and 
accurately  obtained  in  order  to  generate  pointing  data  for  subsequent 
tracking  stations  and  to  determine  contingency  actions  if  the  orbit 
is  unsatisfactory. 
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A  recent  development  in  estimation  theory  is  the  application  of 
the  Extended  Kalman  filter  to  the  early-orbit  determination  problem 
(Ref  3).  The  Extended  Kalman  filter  yields  estimates  of  the  orbit 
states  and  parameters  via  sequential  processing  of  a  set  of  noisy 
observations . 

The  purpose  of  this  study  is  to  compare  the  accuracy  and  ease  of 
application  of  the  Weighted  Least  Squares  filter  and  the  Extended 
Kalman  filter  when  applied  to  the  early-orbit  determination  problem. 

The  application  of  the  Weighted  Least  Squares  filter  to  orbit  deter¬ 
mination  is  described  in  Reference  7.  Previous  studies  considered  the 
application  of  the  Extended  Kalman  filter  to  single  orbit  (Ref  3)  and 
multiple  orbit  (Ref  8)  radar  tracking  observations  of  actual  earth 
satellites.  The  radar  tracking  data  for  these  Kalman  filter  studies 
were  obtained  from  the  Space  Detection  and  Tracking  System  (5PADATS). 
Although  the  results  of  these  two  studies  using  the  Kalman  filter 
compared  favorably  with  the  orbit  parameters  determined  by  SPADATS 
using  a  Weighted  Least  Squares  filter,  the  results  did  not  give  an 
absolute  comparison  of  the  Kalman  filter  with  the  WLS  filter  as  the 
true  orbit  parameters  were  unknown.  This  uncertainty  in  orbit 
parameters  is  due  to  the  fact  that  an  actual  earth  satellite  is 
subjected  to  perturbations  such  as  the  effect  of  a  non-spherical 
earth,  atmospheric  drag  and  magnetic  field  drag.  These  small  perturba¬ 
tions  and  small  long  term  changes  in  orbit  cause  changes  which  were 
not  accurately  modeled  by  the  Extended  Kalman  filter  nor  by  the  WLS 
filter. 

Therefore,  in  this  study,  a  true  two-body  Kcplerian  orbit  is 
simulated,  observations  are  computed  and  contaminated  with  noise,  and 
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the  orbit  parameters  estimated  by  each  filler  are  compared.  In  this 
way,  both  an  absolute  and  relative  coinpn:  i  oi.  can  be  made  of  the 
actual  orbit  parameters  and  estimated  orbit  parameters  produced  by  the 
two  filters.  Hero  the  WLS  filter  as  described  in  Reference  7  and  the 
Extended  Kalman  filter  as  described  in  Reference  4  are  used. 

The  following  assumptions  are  made  concerning  the  two-body 
Keplerian  orbit: 

1.  The  satellite  is  a  point  mans  under  the  gravitational 
attraction  of  a  spherical  rotating  earth  with  negligible  atmospheric 
drag  and  magnetic  field  drag. 

2.  The  satellite  is  non- thrusting. 

The  following  assumptions  are  made  concerning  the  observations  of 
the  orbit: 

1.  Uncertainties  in  latitude,  longitude  and  height  of  the 
tracking  sensors  are  assumed  negligible. 

2.  The  noise  in  the  observation  is  assumed  to  be  an  additive, 
zero-mean,  white,  Gaussian  sequence. 

This  report  is  divided  into  seven  chapters  plus  throe  appendices. 
Chapter  II  includes  an  explanation  of  the  coordinate  systems  used,  the 
required  coordinate  transformations,  a  description  of  the  equations  of 
motion,  the  method  of  determining  the  .initial  conditions,  and  the 
method  used  to  generate  the  noisy  tracking  data.  In  Chapter  III,  the 
Extended  Kalman  filter  equations  and  the  linearization  procedure  used 
by  this  filter  are  described.  Chapter  IV  describes  the  Weighted  Least 
Squares  method  and  the  linearization  procedure  used  in  this  filter. 
Chapter  V  gives  a  comparison  of  the  application  of  the  two  filters  and 
details  their  relative  advantages  and  disadvantages.  A  description  of 
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the  data  and  an  analysis  of  the  results  obtained  are  in  Chapter  VI 
followed  Ly  tho  conclusions  and  recommendations  in  Chapter  VII.  The 
equations  for  determining  tho  initial  position  and  velocity  are 
described  in  Appendix  A;  Appendix  B  contains  a  description  of  the 
system  sensitivity  matrix  and  the  measurement  matrix  description  is 
presented  in  Appendix  C. 


4 


GA/LE/72S-1 


1 1 .  The  Model 

Coordinate  Systems 

Three  coordinate  systems  are  involved  in  the  orbit  determination 
problem:  the  inertial  geocentric  system  and  the  rotating  geocentric 

system  have  as  their  origin  the  center  of  the  earth,  while  the  rotating 
topocentric  system  has  its  origin  on  the  surface  of  the  earth  at  the 
location  of  a  tracking  station.  The  three  systems  are  shown  in 
Figure  1  and  are  described  in  detail  below. 

Inertial  Coordinate  System.  The  center  of  the  inertial  coordinate 
system  is  the  center  of  the  earth.  The  principal  axis,  Xj,  is  directed 
toward  the  first  point  of  Aries,  y.  The  Yj  axis  is  directed  perpen¬ 
dicular  to  Xj  in  the  equatorial  plane,  while  Zj  is  toward  the  earth's 
north  pole.  The  equations  of  motion  employed  in  the  VILS  method  are 
written  in  this  system  because  orbital  motion  is  most  readily  expressed 
in  inertial  space. 

rotating  Geocentric  System.  As  with  the  inertial  system,  the 
center  of  the  earth  is  the  origin  of  the  rotating  geocentric  system. 

The  principal  axis,  Xg,  is  toward  the  earth's  prime  meridian,  Greenwich, 
and  Yg  and  Zg  are  defined  similarly  to  Yj  and  Zj.  The  difference 
between  the  two  systems  is  that  the  rotating  system  rotates  about  the 
Zg  axis  at  earth  rotation,  ij.  Thu  equations  of  motion  used  in  the 
Kalman  filter  are  expressed  in  this  system  in  order’  to  simplify  the 
transformation  to  the  topocentric  system. 

Rotating  Topocentric  Gy  ; t c m .  The  origin  of  the  rotating  topo¬ 
centric  system  is  the  antenna  of  the  tracking  station.  The  principal 
axis,  X-p,  is  toward  local  North  in  the  plane  tangent  to  the  earth  at 


5 


GA/EE/72S-1 


c 


Z£  6  ZR 


Figure  1.  Earth  Centered  Inertial  and  Rotating  and  Station 
Centered  Topocentric  Coordinate  Systems 
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the  station.  The  Z'g  axis  is  straight  up  from  the  station,  while  the 
Y-p  axis  is  toward  local  West  in  the  tangent  plane.  The  actual  and 
computed  values  of  the  measurements  are  expressed  in  this  frame. 


Coordinate  Transformations 

Transformations  between  any  two  reference  systems  are  readily 
made  with  two  matrices:  the  transformation  matrix  between  the  inertial 
and  rotating  geocentric  frames,  and  the  transformation  matrix  between 
the  rotating  geocentric  and  topocentric  frames.  The  matrices  are  given 
below  with  the  angles  as  shown  in  Figure  1. 


xT 

-sin($)cos(  0)  -sin( '>)sin(  0)  cos  ( <> ) 

:<R 

yT 

= 

sin(G)  -cos(0)  0 

yr 

ZT 
.  — 

cos(<fi)coj(G)  cos(  )  s  in  (  0)  sin(b) 

ZR 

XR 

cos(a)t) 

sin(wt)  0 

XI 

vr 

= 

-sin(tot) 

cos  ( (lit )  0 

yi 

ZR 

0 

0  1 

zi 

_  -4 

(2) 


Since  these  are  orthogonal  transformations,  any  other  necessary 
relationships  may  bo  determined  by  multiplying  and/or  transposing  these 
two  matrices. 


Equations  of  Motion 

The  equations  of  motion  of  a  satellite  about  a  spherical  earth 
are  a  set  of  three  second  order  non- linear  differential  equations.  In 
inertial  coordinates  they  arc 
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while  in  earth-centered  rotating  coordinates  they  are 


-yXR  5 

*r  -  +  2t°yR +  w  XR 


yR  =  -  2wxr  +  w2yR 


These  second  order  equations  can  be  reduced  to  first  order 
equations  by  the  use  of  a  state  vector  formulation.  The  state  vector, 

x,  is  defined  as 


Xj  =  x 


Xt*  =  X 


x2  =  y 


x5  =  y 


X(i  =  7. 


The  equations 
equation 


of  motion  then  become  the  6x1  vector  differential 


x(t)  =  f(x(t)> 
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where  f  is  a  6x1  vector  valued  function.  In  inertial  coordinates  the 
state  differential  equation  is 


Xj(  t)  =  fj(Xj(  t)  )  = 


Xu 


x5 


x6 


-UXi 


(11) 


and  in  rotating  coordinates  the  state  differential  equation  is 


Xu 

X5 


*R<t>  =  fR(xK^>> 


x& 

ri^i  + 

r3 

2wx;>  + 

-uy-2 

r3 

2uxi,  +  w2x  2 

-11x3 

Tr 

U 


(12) 
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Initial  Conditions 

In  applying  KLS  and  the  Kalman  filter  to  the  orbit  problem,  an 
initial  or  nominal  state  vector  must  be  available.  Since  the  early- 
orbit  determination  problem  is  very  important  in  non-nominal  launches, 
it  was  decided  to  determine  the  initial  state  from  the  first  set  of 
input  data  rather  than  use  a  nominal  vector.  This  data  set  includes 
range,  range  rate,  azimuth,  elevation,  azimuth  rate,  and  elevation 
rate.  The  angular  rates  are  needed  to  complete  the  set  of  six 
independent  constants  for  the  equations  of  motion.  The  initial 
position  and  velocity  vectors  in  the  rotating  coordinate  system  are 
then  determined  by  the  vector  equations  below  and  the  exact  method 
(Ref  3:03-30)  is  presented  in  Appendix  A. 

r  -  pL  +  R  (13) 

r  =  pb_  pb  (14) 

This  initial  state  vector  is  then  transformed  to  an  inertial  vector 
for  use  in  the  WLS  program. 

Data  Generation 

To  generate  noisy  data  for  the  two  orbit  determination  schemes, 
the  inertial  equations  of  motion  are  integrated  through  the  desired 
time  interval  with  the  starting  position  and  velocity  chosen  from  a 
desired  nominal  orbit.  As  the  satellite  passes  over  the  hypothetical 
station  location,  nominal  tracking  data  is  produced  via  a  coordinate 
transformation  from  .inertial  to  rotating  geocentric  and  topoccntric 
coordinates.  Once  x^,  y^,  zR>  x^,  y^,  Zg,  x-p  yp  and  z-p  are  known, 
the  nominal  observation  data  are  determined  by 
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P  =  C(xR-xs)2  t  (yK-ys)2  +  (xp-:’.s)2lI/2  (lb) 
p  =  “  C(xR-xs)xR  +  (yR-yS)yR  +  <zr-zs)zr]  (16) 
a  =  tan"1  [-yT/x't<]  (17) 
e  =  tan"1  [zT/(xp  +  y2)^2]  (18) 


where  Xg,  yg,  and  ?.^  are  the  station  coordinates  in  the  rotating 
geocentric  system.  If  the  elevation  is  above  a  threshold  value,  the 
observation  data  are  stored  along  with  tiir.e  Iren  to*  When  the 
integration  has  reached  t{ ,  the  table  of  observations  is  contaminated 
by  zero-mean  Gaussian  noise. 

The  random  number  generator  on  the  CDC  0600  produces  a  uniform 
random  number  sequence  between  0  and  1,  thus  the  mean  of  this  sequence 
is  0.5  and  the  variance  is  1/12.  To  convert  this  sequence  to  an 
approximate  normal  sequence  of  mean  zero  and  variance  one,  the  Central 
Limit  Theorem  is  used.  This  theorem  applied  to  a  uniform  sequence 
states  that 


N  ( 0 , 1 ) 


(10) 


where  Sn  is  the  sum  of  n  uniform  numbers,  p  is  the  uniform  mean,  o  is 
the  uniform  standard  deviation,  and  is  the  normal  distribution 
function.  In  n=12,  the  equation  reduces  to 
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12 

J[  x.  -  6  «  11(0,1)  (20) 

i  =  l 

and  a  normal  die tribution  curve  is  approximated  by  a  12th  order 
polynomial.  This  is  the  equation  used  to  produce  the  approximate 
zero-mean  Gaussian  sequence.  To  change  the  variance,  each  random 
number  is  multiplied  by  the  desired  standard  deviation,  producing  a 
i! ( 0 , o 2 )  distribution.  If  started  at  the  same  place,  this  random 
number  generator  will  always  generate  the  same  sequence  realization; 
thus,  to  obtain  different  realizations  the  generator  is  started  with 
a  random  initial  condition.  This  normal  distribution  approximation 
suffers  from  one  disadvantage  in  that  only  a  range  is  produced, 
thus  the  "tails"  of  the  curve  are  not  generated.  The  impact  of  this 
on  the  orbit  problem  is  that  very  large  data  errors  are  not  encountered 
and  the  filters  may  periorm  better  than  they  would  with  normally 


distributed  data. 
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III.  The  Kalman  F i Ltcr  Equa t tops 

Measurement  Model 

The  Kalman  filter  as  described  in  Reference  4  processes  noisy 
linear  data  in  a  recursive  manner  to  produce  a  linear  unbiased  minimum 
variance  estimate  of  the  state  of  a  linear  dynamic  system.  The  data 
consist  of  a  linear  function  of  the  state  plus  additive  Gaussian 
noise . 

In  the  orbit  determination  problem,  the  data  are  ground  radar 
measurements  range  (p),  range  rate  (p),  azimuth  (a),  and  elevation  (e) 
which,  as  seen  in  equations  (15-18)  are  not  linearly  related  to  the 
state  vector  and  the  station  position  vector.  The  measurement  model 
may  be  expressed  as  a  discrete  non-linear  vector  equation 

z(k)  =  h[x(k)j  +  v(k)  (21) 

where  x  ia  a  4x1  vector  consisting  of  p,  p,  a,  and  o;  h  is  a  non-linear 
vector  valued  function;  and  v  is  a  v/hito  Gaussian  sequence.  It  is 
assumed  that  the  noise  has  mean 

E[v(k)]  =  0  (22) 

in  this  case  but  this  assumption  is  not  necessary  in  the  general 
problem.  The  noise  covariance  .is  assumed  to  be 

E[v(k)vT(k)]  =  R(k)  6-j  (23) 

where  R(k)  is  a  known  positive  semi-definite  matrix. 


.1.3 
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Filter  Equations 

Since  the  Kalman  filter  assumes  a  linear  system  ind  a  linear 
measurement  model,  the  filter  cannot  be  applied  directly  to  estimating 
the  state  of  a  satellite  orbit  based  on  a  sequence  of  noisy  radar 
observations.  Equations  (12)  and  (21)  are  the  non-linear  state  and 
measurement  equations  and  must  be  approximated  by  a  linearised  set  of 
equations  before  the  linear  filter  can  be  used.  The  two  sets  of 
equations  are  expanded  in  a  Taylor  series  about  a  suitable  point  and 
first-order  terms  are  retained.  The  linearized  equations  are 

x(t)  et  Fx(  t)  (24) 


and 


z(k)  =  Mx(k)  +  v(k) 


where  F  is  the  6x6  State  Sensitivity  Matrix 
3f(x) 


F(x)  = 


3x 


x  =  x(k+l|k) 


and  M  is  the  4x6  Measurement  Matrix 


(25) 


(26) 


M(x) 


3h(x) 

3x 


x  =  k(k+.l.|k) 


(27) 


The  exact  form  of  these  matrices  is  developed  in  Appendices  13  and  C. 

The  nominal  point  about  which  the  non-linear  equations  are 
expanded  is  the  current  state  estimate  x  based  on  the  previous  measure¬ 
ment  at  time  k  and  integrated  forward  via  equation  (12)  to  time  k+1. 

As  a  result  of  updating  the  nominal  at  each  time  point,  initial  errors 


14 


GA/EE/72S-1 


are  not  allowed  to  propagate  through  time  and  the  linearity  assumptions 
should  be  more  valid  as  the  filter  progresses  (Ref  4:276). 

Now  that  the  state  and  measurement  equations  are  in  a  linear  form, 
the  recursive  Kalman  filter  equations  can  be  applied.  The  filter 
consists  of  a  predictor  in  which  the  estimated  state  and  state  error 
covariance  are  integrated  forward  to  the  next  observation  time  point 
and  a  corrector  in  which  the  state  estimate  and  estimated  state  error 
covariance  are  corrected  by  the  new  observation  data. 

Prediction.  Prediction  from  time  t]<  to  time  t^+1  is  accomplished 
by  numerically  integrating  the  state  vector  differential  equation  (12). 
In  discrete  form  this  equation  becomes 


a(k+l|k)  =  ft(k 


f[x(k |k)]dt 


(28) 


The  state  error  covariance  at  time  tj<+^  based  on  data  at  time  t^  is 

P(k+l|k)  =  >!>(k+l,k)P(k|kHT(k+l,k)  (29) 

where  <I>(k+l,k)  is  the  State  Transition  Matrix: 

as(k+i) 

<J>(k+l,k)  =  — -  (30) 

3x(k) 

<!>(k+l,k)  is  determined  by  integrating  the  linear  differential  matrix 
equation 

I 

(31) 


•j>(k+l,k)  -  F[ft(k|k)]<:(ktl,k) 
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The  error  covariance  of  the  corrected  state  estimate,  R(k+l|k+l),  is 
determined  by  the  equation 

P(k+1  |k+l)  =  Ll-K(ktl)MjP(kH|k)Ll-K(k+i)M]T  + 

K(k+l)RK(k+l)T  (39) 

where  M  is  again  calculated  at  time  t^^. 

To  start  the  process,  values  of  jUo|o),  P(o|o),  and  K  must  be 
available.  The  initial  state  estimate,  k_(o|o),  is  computed  by  the 
method  in  Chapter  II;  the  initial  state  error  covariance  matrix, 

P(0 | 0) ,  must  be  estimated  a  priori;  and  the  measurement  error  covari¬ 
ance  matrix,  R,  is  known  from  antenna  tests,  or,  when  the  input  data 
are  simulated,  from  the  variances  of  the  noise  added  to  the  measurements. 

The  filter  used  in  this  manner  is  referred  to  as  the  Extended 
Kalman  filter  due  to  the  fact  that  the  equations  are  linearized  about 
the  current  state  estimate,  £,(k+l|k).  The  prediction  process  does  not 
involve  any  linear  assumptions,  as  the  non-linear  equations  (12)  and 
(33)  are  used  to  predict  the  state  and  the  measurement  and  the  state 
error  covariance  is  predicted  by  the  linear  matrix  differential 
equation  (29).  The  linearization  procedure  is  only  used  in  the 
calculation  of  the  corrected  state  estimate  and  corrected  state  error 


17 


GA/EE/72S-1 

covariance,  i.e.,  linearization  involves  the  state  errors  rather  than 
the  states  themselves. 

A  derivation  of  the  Extended  Kalman  filler  is  found  in  Reference  4 
and  a  flow  diagram  for  the  filter  is  presented  in  Figure  2. 
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IV.  Weighted  Leant  Squares 


Least  Squares  Theory 

The  goal  of  the  method  of  Weighted  Leant  Squares  when  applied  to 
orbit  determination  is  to  determine  the  values  of  the  six  orbit 
parameters  at  some  epoch  time,  tg,  given  sets  of  observation  data  such 
that  the  cost  function 

J  =  III  -  £(£o)|lw  3  Cz  -  £<£o>JT  WCZ  “  (,4°) 

is  a  minimum.  W  is  an  appropriate  weighting  matrix;  y  is  the  observa¬ 
tion  vector  consisting  of  range,  range  rate,  azimuth,  and  elevation  at 
each  discrete  time  tfc,  (tg  <  tj^  5  tf),  and  £  is  the  predicted  observa¬ 
tion  vector  based  on  the  current  state  vector  estimate  £  at  tg.  Thus 
y  is  the  hn*l  measurement  vector 


— t0 

Sxi 

-t2 


and  £  is  the  hn*l  estimated  measurement  vector 
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(42) 


'31 


The  value  of  X£  which  minimizes  the  cost  function  is  the  "best" 
estimate  in  the  least  squares  sense. 

Since  the  predicted  observation  vector  £  is  a  non-lincar  function 
of  the  states,  the  equation  relating  them  must  be  linearized  by 
expanding  £  in  a  Taylor  series  about  and  dropping  higher  order  terms. 
The  cost  function  becomes 


J  =  I  lAv|  | J  =  |  lA^  -  AAxqHJ 


(43) 


where  A  is  the  4n*6  Measurement  Sensitivity  Matrix 


(44) 


and 


=  y.  -  £(£o) 


(45) 
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Differentiating  equation  (43)  v;ith  respect  to  the  state  correction 
vector,  Axq,  and  setting  the  result  equal  to  zero: 

2[ATWAAx0  -  A'V]T  =0  (46) 

This  equation  is  then  solved  for  Axq  ,  yielding  the  normal  matrix 
equation  for  the  WLS  problem: 

AtWAAx0  =  ATWAy_  •  (47) 

or 

Ax0  -  [ATWA]'1  ATWAZ  (48) 

The  differences,  Ay_,  are  the  residuals  of  the  measurements  and  are  the 
result  of  the  random  observation  noise ,  incorrect  values  of  the 
estimated  state  q  ,  anQ  computation  errors.  Since  the  original  WLS 
problem  has  been  reduced  to  the  linearized  problem  of  determining  a 
correction  vector  A><o  such  that 

IUz-AAxoHJ  (49) 

is  minimized,  the  solution  is  only  approximate.  Therefore,  an  iterative 
process  must  be  used  until  the  residuals  are  reduced  to  some  minimum 
(Ref  7:2-4). 

Convergence  Criteria 

In  order  to  determine  when  to  stop  the  iteration  process,  some 
convergence  criteria  must  be  employed.  Using  the  corrected  state 
vector  to  predict  a  now  set  of  weighted  residuals 
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I  lAXP|  lw  =  I  lAX  "  AA*ol  lj  (5°) 

The  expression  | |Ay^| is  then  the  RMS  of  the  predicted  residuals. 

The  least  squares  process  is  said  to  converge  if  either 

IMIw-  Ill’ll* 

-  <  ej  (51) 

IWI„ 

or 

I  I  All  lw  -  I  lAy'U  <  e2  (52) 

where  |  |  Ay  |  |  is  the  RMS  of  the  current  weighted  residuals  and  ej  and 
e2  are  some  small  constants  (Ref  7:20).  The  values  chosen  for  ei  and 
e2  for  this  study  were  both  2xlQ“'f. 

Sensitivity  Matrix 

The  Measurement  Sensitivity  Matrix,  A,  may  be  calculated  using 
the  chain  rule: 


3xt 

Dxq  3_x t 


(53) 


and  defining 


(54) 
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Statistical  Aspects 

Until  now,  nothing  has  been  said  about  either  statistical  theory 
or  the  contents  of  the  weighting  matrix,  W.  For  many  years,  this 
matrix  was  guessed  for  each  problem  and  was  used  as  a  "fudge  factor" 
since  it  seemed  reasonable  to  weight  "reliable"  data  more  heavily 
than  "unreliable"  data.  However,  recent  statistical  theory  shows  that 
if  this  weighting  matrix  is  chosen  to  be  the  inverse  of  the  observation 
error  covariance,  R,  and  the  observation  errors  have  a  r.ero-mean,  then 
the  matrix 

P  =  ATWA  --  AtR-1A  (60) 

is  the  state  error  covariance  matrix.  In  addition,  the  estimated 
state  error,  Axq ,  which  results  from  equation  (48)  using  the  covari¬ 
ance  matrix,  P,  is  the  linear  unbiased  minimum  variance  estimate  of 
the  state  error.  If  the  observation  errors  have  a  joint  Gaussian 
distribution,  then  the  state  error  estimate  is  also  the  maximum 
likelihood  estimate  (Ref  7:25-26). 

Computational  Procedure 

The  procedure  used  to  determine  the  satellite's  orbit  parameters 
using  WLS  are  as  follows: 

1.  Set  A  =  $.  Input  R  and  £q  where  vq  consists  of  p,  p,  a, 
and  e. 

2.  Determine  an  initial  nominal  state  vector,  ,  from  the  first 
data  set  (p,  p,  a,  e,  a,  and  e). 

3.  Integrate  the  trajectory  and  the  State  Transition  Matrix  from 
to  to  tf. 
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(n)  At  each  time  point,  tfc ,  calculate  and  store 


h  =  £<$o) 

the  predicted  measurement 

Mfc,  the  Measurement  Matrix 

Ak  =  Ak-.l  +  Mk*<tk»t0> 


A^k  r-  -  i<io>k 


the  measurement  residual. 
4.  When  t  =  tf,  calculate 
P  =  ATR_iA 


the  state  error  covariance,  and 


llAzl  lu 


the  weighted  residual  vector  RMS. 


5.  Invert  the  covariance  matrix  and  calculate 


Ax0  =  P“1AtR“1A£ 


then 


£0  =  £0  +  AiiO 


(61) 


(62) 

(63) 


(64) 


(65) 


(66) 


(67) 


(68) 
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6.  Calculate  the  predicted  RMS 

llAy.FMw  =  AAx0|  |w31/2  (69) 

and  determine  whether  or  not  the  convergence  criteria  have  been  met. 

If  they  have  not,  return  to  step  3  with  the  new  Rr, .  If  the  convergence 
criteria  have  been  met,  stop.  If  the  predicted  RMS  is  greater  than 
the  current  RMS,  or  if  the  current  RMS  is  greater  than  the  old  RMS, 
the  process  is  diverging,  indicating  some  problem  with  either  the  data 
or  the  system  model.  In  this  case,  the  process  stops  and  all  of  the 
data  must  be  examined.  One  type  of  data  which  can  cause  divergence 
is  an  "outlier"-~a  data  point  very  far  from  the  other  points.  This 
point  can  be  spotted  by  looking  at  the  residuals  of  the  measurements 
and  must  be  manually  thrown  out.  The  WLS  procedure  would  then  be 
restarted  without  the  outliers. 

A  computer  flowchart  for  the  WLS  method  is  shown  in  Figure  3. 
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V .  Comparison  of  the  Opera t ion  of  th  Two  Pi  1  tors 
Theory  of  the  Two  Methods 

It  was  stated  in  Chapters  III  and  IV  and  shown  in  References  4  and 
7  that  both  of  the  methods  produce  a  linear  unbiased  minimum  variance 
estimate  of  the  state  of  a  linear  system  given  linear  measurements  with 
additive  noise.  Reference  7  (page  157)  lists  a  number  of  authors  who 
have  shown  that  when  the  Weighted  Least  Squares  method  is  applied  in  a 
recursive  manner,  the  resulting  equations  are  exactly  the  same  as  those 
in  the  Kalman  filter.  Therefore,  the  only  difference  between  the  two 
methods  compared  in  this  paper  is  that  one  is  a  batch  processor  while 
the  other  is  recursive.  Doth  are  linear  processors  and  both  linearize 
the  earrors  in  the  orbit  problem.  The  extended  Kalman  filter  linearizes 
about  the  current  state  estimate  at  time  t;<}  while  WLS  linearizes 
about  the  state  estimate  at  Iq. 

Advantages  and  Disadvantages 

Since  the  two  methods  arc  of  the  same  form,  the  only  theoretical 
comparisons  which  can  be  made  are  those  between  batch  and  recursive 
processing.  The  following  are  a  few  of  the  theoretical  advantages  and 
disadvantages  of  the  Hxtended  KaJrnnn  filter  as  compared  with  the 
method  of  Weighted  Least  Squares. 

A.  Since  the  filter  linearizes  about  the  current  state,  a  large 
initial  error  plus  many  observations  over  a  long  time  span  would 
probably  not  cause  the  filter  to  diverge  while  WLS  may  diverge.  This 
is  because  WLS  would  use  the  erroneous  initial  orbit  to  predict  for 
the  entire  span  and  errors  at  later  times  may  be  so  large  that  the 
linearity  assumptions  are  violated. 
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B.  With  a  largo  amount  ot  data,  both  methods  tend  to  ignore 
changes  In  new  data.  In  WLS,  the  new  data  are  overcome  by  the 
preponderance  of  older  data  and  in  the  filter,  once  it  "learns"  the 
model,  the  Kalman  gain  becomes  so  small  that  new  data  are  under¬ 
weighted.  The  solution  in  both  canes  is  to  restart  the  process  and 
leave  out  most  of  the  old  data. 

C.  The  filter  data  storage  is  small  while  the  WLS  residual 
vector  must  be  stored  and  may  be  very  large. 

D.  The  filter  processing  requirements  are  quite  a  bit  Larger 
than  the  requirements  of  WLS. 

E.  An  stated  in  Reference  1  (page  390),  the  filter  processor  is 
slower  than  one  iteration  of  WLS,  however,  WLS  usually  requires  more 
than  one  iteration,  so  they  both  require  similar  central  processor 
time. 

F.  The  filter  produces  an  estimate  of  the  state,  the  state  error, 
and  the  state  error  covariance  at  each  point  in  the  process.  The 
estimated  error  and  error  covariance  are  extremely  useful  in  analysis 
of  the  process  as  any  outliers  in  the  data  are  very  apparent,  so  that 
real-time  control  of  the  data  is  possible.  The  WLS  method  gives  a 
weighted  RMS  and  predictei  RMS  in  addition  to  the  state  error  covari¬ 
ance,  but  these  only  givj  an  indication  of  what  the  processor  did  after 
an  iteration.  The  whole  residual  vector  must  be  examined  to  discover 
outliers,  and,  if  the  process  diverged,  there  is  usually  no  clear 
reason  for  the  problem,  nor  is  there  shown  any  time  point  at  which  the 
process  started  to  blow  up. 
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c 


This  time  history  output  of  the  filter  is  probably  the  most 
important  advantage  compared  to  *LS.  It  can  Le  used  to  determine  the 
exact  time  at  which  a  problem  occurs,  and  if  the  process  starts  to 
diverge  at  any  point,  it  may  be  restarted  after  that  point  with  only 
the  initial  guess  of  P(k|k)  and  £(k)  needed. 
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VI.  Rc-sM.lt:; 

In  order  to  compare  the  performance  ot  the  two  filters  in  the 
early-orbit  determination  problem,  noisy  observations  produced  by  the 
generation  program  and  necessary  station  data  die  input  t:o  the  two 
filtering  programs.  The  comparison  then  consists  of  two  parts.  First, 
to  compare  the  application  of  the  two  methods  to  this  problem,  a  single 
pass  of  observation  data  is  processed  by  each  program.  The  same  data 
set  is  input  to  each  computer  program  so  that  the  results  can  be 
compared.  The  resulting  outputs  of  the  two  filters  are  than  compared 
to  determine  their  suitability  to  the  orbit  determination  problem  and 
to  show  the  type  of  information  they  each  provide.  The  titters  are 
rerun  for  varying  numbers  of  observations  so  that  the  effects  of 
observation  truncation  may  be  seen. 

Next,  to  obtain  a  comparison  of  the  two  filters'  absolute 
accuracies,  their  radial  and  velocity  estimates  for  various  numbers  of 
observations  are  compared  to  the  actual  orbit  radii  and  velocities. 

This  provides  an  absolute  accuracy  comparison  of  the  two  methods. 

If  only  one  sot  of  observation.',  were  used  in  making  this  absolute 
comparison,  the  results  would  not  be  valid.  This  is  because  the 
observation  sequence  is  finite  and  the  additive  noise  is  a  truncation 
of  an  iniinite  noise  sequ  nee  which  ha;  Gaussian  itatistics.  Therefore, 
the  statistics  of  any  single  truncated  real! sal  ion  of  the  .infinite 
sequence  will  probably  not  have  the  desired  Gaussian  statistics.  A 
reasonable  approximation  to  the  Gaussian  statistics  on  the  additive 
noise  is  obtained  by  taking  an  ensemble  average  of  the  filters'  outputs 
for  50  different  truncated  r  indent  sequences  u  ed  fox’  the  additive 
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noi so.  Th*j  results  to  be  compared  are  then  the  mean  I'adial  and 
velocity  errors  of  the  50  runs.  Figure  4  shows  a  flow  diagram  of  the 
accuracy  comparison  process.  Fifty  data  sots  of  58  observations  each 
are  produced,  radial  and  velocity  vectors  are  estimated  by  the  two 
filters  and  the  magnitudes  of  those  vectors  are  compared  with  their 
actual  values  from  the  nominal  orbit. 

Data 

The  actual  orbit  parameters  of  the  two-body  orbit  are  given  in 
Table  I  along  with  the  location  of  the  hypothetical  tracking  station 
and  the  sigmas  of  the  noisy  observations.  A  medium-altitude  orbit 
was  chosen  because  atmospheric  drag  and  solar  pressure  are  small  in 
this  region  and  the  assumptions  of  a  two-body  problem  with  no  drag 
would  be  mere  realistic. 

Since  the  problem  is  an  early-orbit  problem,  only  one  tracking 
pass  was  used.  This  particular  data  set  simulates  a  near-polar 
satellite  launched  at  Vandenbcrg  Air  Torcc  Base  passing  over  Shemya 
Tracking  Station  on  the  first  orbit. 

Filter  Output  Comparison 

Table  II  shows  the  error  magnitudes  predicted  by  the  Kalman 
filter,  the  standard  deviations  predicted  by  each  of  the  methods  and 
the  actual  radial  and  velocity  error  magnitudes  for  a  various  number 
of  observations.  Since  these  estimated  value  s  are  for  only  one  pass 
of  data  which  ranged  from  4  to  50  observations ,  the  actual  numbers 
give  an  indication  of  the  order  of  magnitude  ol  the  outputs  of  the 
two  programs.  The  main  result  of  this  comparison  lies  in  the  trends 
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Tdblo  I 

Actual  Orbit  and  Sensor  Data 


True  Two-Budy  Kcplcrian  Orbit  I 

Parameters 

Period 

Inclination 

Eccentrici ty 

06.38205  Min 

56.0713  Deg 

0.00312689 

Longitude  of 

Argument  of 

Perigee 

Ascending  tlode 

Perigee 

Height 

203.3325  Deg 

21  8. 1018  Deg 

306.2086  NM 

Radius  at  tQ 

Velocity  at  Lq 

6003973.2  !i 

75113.61  M/S 

Tracking  Station  Sensor  Data 

Latitude 

Longitude 

Height 

52.73267  N 

174.1023  E 

0.0 

Azimuth 

Sigma 

(Deg) 

Elevation 

Sigrr.a 

(Dog) 

Range 

Sigma 

(M) 

Range  Rate 
Sigma 
(M/S) 

0.02 

0.02 

100.0 

1.0 

Least  Squares  Radial  and  Velocity  Errors  -  Single  : 
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of  the  error:;  and  standard  deviations  as  the  number  of  observations 
is  varied. 

As  seen  in  Table  II,  the  estimated  sigmas  produced  by  each  method 
are  fairly  close  but  the  Least  Squares  radial  and  velocity  sigmas 
decrease  monotonically ,  while  the  Kalman  filter's  radial  sigma  reaches 
a  low  at  around  30  tracking  points  and  then  starts  climbing  slowly. 
This  is  probably  due  to  the  accumulation  of  computer  round-off  error. 
Since  the  Kalman  filter  computes  the  covariance  matrix  in  a  recursive 
manner,  any  errors  will  bo  additive  and  will  finally  affect  the 
significant  digits  of  the  estimates  covariance  and,  subsequently,  of 
the  estimated  state.  Since  the  radial  variance  is  so  much  larger  than 
the  velocity  variance,  this  effect  would  bo  observed  first  in  the 
radial  sigma,  as  seen  in  Table  II.  The  WLS  process  computes  the 
estimated  state  and  state  covariance  independently  from  one  iteration 
to  the  next,  so  round-off  errors  cannot  propagate  from  one  calculation 
to  the  next.  Thus  the  radial  and  velocity  sigmas  produced  by  this 
filter  do  not  show  any  growth. 

The  estimated  state  error  produced  by  the  filter  has  no  counter¬ 
part  in  WLS.  Since  WLS  does  all  of  its  state  calculations  at  to,  no 
time  history  of  state  errors  is  generated.  As  seen  in  Table  II,  these 
estimated  error  values  are  not  close  to  the  actual  errors  in  most 
cases  but  they  do  give  an  indication  of  the  order  of  magnitude  of  the 
actual  errors  and  of  the  stability  of  the  process.  Whenever  the 
estimated  error  is  greater  than  the  estimated  sigma  a  majority  of  the 
time  the  filter  is  diverging.  When  these  values  are  plotted,  as  in 
Figures  5  and  6,  they  give  a  good  picture  of  the  operation  of  the 
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filter.  This  is  the  "big  picture"  which  Wi.b  cannot  j  reduce,  figures 
5  and  G  show  that  the  filter  has  no  tendencies  toward  divergence  in 
this  case.  No  divergent  examples  occurred  in  this  study,  and  none 
were  expected  since  the  model  and  the  actual  system  were  the  same 
non-linear  system  of  equations.  In  addition,  the  noisy  data  used  in 
this  problem  wore  reasonably  well-behaved  since  they  contained  no 
additive  noise  larger  than  six  sigma.  In  a  real  situation,  a  pre¬ 
processor  would  be  used  to  remove  extremely  noisy  data,  so  the  results 
of  this  study  should  be  applicable  to  real  systems  where  the  actual 
model  is  known  fairly  well. 

Orbit  Parameters .  When  an  early-orbit  determination  program  is 
actually  used,  the  program's  output  which  has  physical  significance 
is  not  the  radius  or  velocity  but  the  orbit  parameters.  Thus,  to 
compare  the  two  filters'  application  to  this  problem  requires  a 
comparison  of  the  estimated  orbit  parameters  with  the  actual  parameters. 

Table  III  shows  the  orbit  parameters  computed  from  the  estimated 
radius  and  velocity  produced  by  each  method  from  various  numbers  of 
observations.  The  method  of  calculating  the  two-body  orbit  parameters 
from  radius  and  velocity  is  given  in  Reference  2  (pages  98-107).  Even 
with  as  few  as  ten  observations,  the  predicted  orbit  parameters  are 
fairly  close  to  the  actual  values,  and  with  the  whole  pass  the  errors 
are  extremely  small. 

The  most  important  element  to  be  computed  in  the  early-orbit 
situation  is  the  period.  This  is  because  the  pointing  data  at  the  next 
station  are  very  sensitive  to  period  errors.  A  period  change  of  only 
10  seconds  could  cause  the  satellite  to  be  outside  the  beam  width  of 
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the  tracking  antenna.  The  period  in  a  simple  function  of  the  radius 
and  velocity,  as  seen  from  the  lollowing  equations. 


period  (seconds) 


(70) 


where 


(71) 


and  k  is  the  Earth  gravitation  constant  in  radians/second.  Thus  the 
sensitivity  of  the  period  to  errors  in  radius  and  velocity  is  a  direct 
function  of  errors  in  each  one.  As  seen  in  Table  III,  the  period 
estimated  by  each  filter  using  58  observations  is  much  less  than  one 
second  from  the  true  value,  hence  the  satellite  should  be  easy  to 
find  at  the  second  and  subsequent  stations. 

Accuracy  Comparison  Using  50  Data  Sets 

Figures  7  and  8  show  the  two  filters'  mean  radial  and  velocity 
error  for  50  runs  versus  the  number  of  observations  while  Figures  9 
and  10  show  the  radial  and  velocity  error  standard  deviations.  The 
errors  are  defined  as  estimated  minus  actual  values.  Figure  9  shows 
that  the  Kalman  filter  has  a  slightly  smaller  radial  error  sigma  when 
less  than  30  observations  arc  processed,  while  Kbd  appeal's  slightly 
bettor  tor  data  amounts  greater  than  that  threshold.  This  appears 
consistent  with  the  theory  since  the  Kalman  filter  "learns"  the  model 
faster  than  the  WLS  filter,  but  once  the  Kalman  filter  errors  read)  a 
minimum,  the  covariance  error  starts  growing  and  affects  the  estimate 
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accuracy.  WLS  does  not  experience  tin.  covariance  error  o  the  actual 
error’s  continue  to  decrease  with  the  number  oi'  observations  processed. 

Figure  10  shows  that  the  velocity  error  distributions  of  both 
methods  are  so  similar  that  they  are  indistinguishable.  Since  the 
velocity  magnitudes  are  so  much  smaller  than  the  radial  magnitudes, 
the  difference  in  tiie  errors  in  velocity  are  not  as  noticeable  as  the 
radial  errors. 

Non-Linearity  Problems.  Figures  7  and  13  reveal  a  problem  with 
the  Kalman  filter.  The  estimates  of  the  two  methods  should  be  unbiased, 
however  the  Kalman  filter's  mean  radial  and  velocity  errors  show  a 
small  negative  bias  after  some  rather  large  initial  errors.  This  bias 
of  approximately  30  meters  in  radius  represents  an  error  of  0.0005% 
and  the  velocity  bias  of  approximately  0.1  m/s  is  an  error  of  0.001%, 
so  it  would  never  be  apparent  when  one  case  is  run  at  a  time.  Only 
when  several  cases  wore  studied  did  the  small  biases  become  apparent. 

The  effect  of  error’s  of  these  magnitudes  is  soon  in  Table  III  in  which 
the  velocity  and  radial  errors  for  50  observations  are  nearly  the  same 
as  the  mean  errors  of  the  50  runs. 

In  the  search  for  some  bias  in  the  computer  program,  it  was  found 
that  the  .initial  velocity  estimated  by  the  initial  condition  module 
was  always  500  to  700  m/s  higher  than  the  actual  velocity.  This  large 
velocity  error  is  nearly  10%  of  the  actual  velocity  and  was  due  to  the 
largo  errors  in  the  angle  and  angular  rate  data.  If  these  errors 
exceed  the  linearity  region  in  the  Kalman  filter  assumptions,  a 
possible  result  is  a  biased  estimate  (Ref  7:343).  To  see  if  this  was 
the  eai.e,  an  initial  velocity  with  error  magnitude  of  approximately 
100  m/s  was  manually  put  into  the  program.  The  resulting  mean  radial 
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and  velocity  error;;  for  the  50  cases  are  shown  in  figui*  11  and  12. 

The  V/LS  mean  errors  were  exactly  the  same  as  in  the  first  50  runs  and 
are  included  so  that  the  change  of  scale  can  be  seen.  Those  figures 
show  that  the  mean  errors  for  the  smaller  amounts  of  data  do  not 
experience  the  large  errors  observed  in  the  original  runs  and  that  the 
negative  bias  observed  has  been  reduced  by  nearly  50%.  The  most 
probable  reason  for  the  large  errors  and  the  bias  is  that  the  linearity 
assumptions  were  violated.  The  Kalman  filter's  sensitivity  to  non- 
linearities  is  discussed  by  Jaxwinski  (Ref  4:348).  This  example  shows 
that  V/LS  is  apparently  not  as  sensitive  to  non-linearity  us  the  Kalman 
filter  since  it  showed  no  biases  or  erratic  results  when  used  with  the 
larger  initial  errors. 

To  determine  non-linear*  effects  in  the  V/LS  filter,  several  rims 
were  made  with  progressively  larger  initial  velocity  errors.  The  WLS 
process  showed  no  increased  error  until  the  initial  velocity  errors 
were  in  the  6000-7000  m/s  range.  Specifically,  the  filter  diverged 
in  all  cases  when  the  errors  reached  6200  m/s  and  then  no  orbital 
elements  were  produced. 

To  see  how  errors  of  this  magnitude  affected  the  Kalman  filter, 
eases  were  run  with  an  initial  velocity  error  of  7500  m/s.  The 
resulting  mean  radial  and  velocity  errors  are  shown  in  Figures  13  and 
14.  For  58  observations,  the  orbital  periods  produced  ranged  between 
86.29536  and  96.30740  minutes.  These  values  reflect  an  error  range 
of  4.53  to  5.19  second.  Thus,  with  initial  error  magnitudes  larger 
than  those  which  cause  the  Weighted  Least  Squares  method  to  diverge 
and  produce  no  usable  .results,  the  Kalman  filter  is  still  able  to  give 
relatively  accurate  estimates  of  the  critical  orbit  element. 
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Computer  HoTiQi'y  and  Time  Rt  quirenents 

The  CDG  660  centr.il  processor  time  inquired  by  each  program  to 
compute  the  final  orbit  wan  similar  in  most  cases,  with  the  Kalman 
filter  being  slightly  faster. 

For  10  data  points,  the  time  per  case  for  the  Kalman  filter  was 
1.02  seconds  and  for  WLS  1.29  seconds.  For  58  data  points,  the  time 
per  case  was  3.0  seconds  for  the  Kalman  filter  and  5.19  seconds  for 
the  Least  Squares  filter. 

The  core  storage  requirements  for  one  pass  of  data  for  the  two 
programs  V'cre  67200  octal  words  for  tl  _>  Kalman  filter  and  60200  octal 
words  for  Least  Squares.  The  majority  of  the  space  used  by  the  Kalman 
filter  was  for  the  filter  itself,  whereas  in  the  WLS  filter  a 
significant  proportion  of  the  memory  was  used  for  residual  and  matrix 
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Figure  7.  Kalman  and  Least  Squares  Mean  Radial  Error  -  50  Runs 
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Figure  8.  Kalman  and  [.east  Squares  Mean  Velocity  Error  -  50  Run 
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Figure  3.  Kalman  and  Least  Squares  Radial  Error 
Standard  Deviation  -  50  Runs 
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Figure  10.  Kalman  and  Least  Squares  Velocity  Error 
Standard  Deviation  -  50  Runs 
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figure  11.  Kalman  and  Least  Squares  Mean  Radial 
Error,  Av0  =  100  m/s 


Figure  12.  Kalman  and  Least  Squares  Mean  Velocity 
Error,  A vq  =  100  m/s 
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VII .  Conclusions  and  Recommendations 

Conclusions 

1.  The  Weighted  Least  Squares  filter  and  the  Extended  Kalman 
filter  produce  extremely  accurate  results  when  applied  to  one-pass 
tracking  data  for  a  medium-altitude  satellite  where  the  model  is  known. 
The  accuracy  of  both  methods  drops  when  fewer  observations  are  used, 
but  with  as  few  as  10  observations  usable  answers  are  still  produced. 

2.  For  the  orbit  determination  problem  considered  the  Extended 
Kalman  filter  is  more  sensitive  to  non-linearities  in  the  model  than 
Weighted  Least  Squares  and  is  more  sensitive  to  initial  condition 
errors.  This  sensitivity  leads  to  relatively  large  mean  errors  in 
the  early  part  of  the  process  and  an  overall  small  bias  as  the  process 
continues.  These  biases  are  not  of  a  magnitude  to  seriously  affect 
the  resulting  orbit  parameters,  but  they  do  grow  to  fairly  large 
magnitudes  when  the  initial  errors  become  very  large.  When  these 
initial  errors  finally  become  large  enough  to  cause  divergence  in  the 
Weighted  Least  Squares  process,  the  Kalman  filter  is  still  able  to 
produce  usable  estimates. 

3.  The  computer  memox’y  and  central  processor  time  requirements 
are  similar  for  both  methods. 

Recommendations 

The  following  topics  are  recommended  for  further  study: 

1.  Apply  the  Extended  Kalman  filter  and  the  Weighted  Least 
Squares  filter,  using  the  same  model,  to  data  with  different  variances 
and  to  data  from  different  orbits  to  continue  the  comparison. 
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2.  Examine  the  present  orbit  model  ujiiv  b'vighted  Least  Squares 
compared  with  the  Extended  Kalman  fiatei  with  smoothing  back  to  tg. 

3.  Compare  the  two  methods  when  applied  to  multiple  passes  of  a 
satellite. 

4.  Examine  the  two  filters  when  applied  to  the  present  model 
with  known  drag  forces  added 

5.  Apply  the  methods  using  the  simple  earth  model  to  actual  data 
and  compare  their  handling  of  model  errors. 

6.  Expand  the  two  systems  to  employ  the  current  known  earth 
model  and  use  them  with  actual  data.  These  results  could  be  compared 
with  results  from  the  previous  recommended  study. 

7.  Combine  the  two  methods  in  a  Limited  Memory  Filter  as 
described  in  Reference  4  (pages  255-250).  This  method  would  combine 
the  best  points  of  each  of  the  two  methods. 
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Appendix  A 

Determination  of  Initial  Conditions 

The  equations  used  to  produce  the  initial  state  estimate  are 
presented  in  this  appendix.  The  initial  condition  module  needs  the 


following  data: 

1.  Tracking  station  height  in  earth  radii:  h 

2.  Tracking  station  latitude  in  degrees:  4> 

3.  Tracking  station  longitude  in  degrees:  0 

4.  Slant  range  in  meters:  p 

5.  Slant  range  rate  in  metei's/second:  p 

6.  Azimuth  from  north,  clockwise,  in  degrees:  a 

7.  Elevation  in  degrees:  e 

8  Azimuth  rate  in  degrees/second:  a 

9.  Elevation  rate  in  dogrees/second:  e 


Station  Coordinates  in  Rotating  Frame 


Xo 

j 

( l+h)cos(  0)cos(  t}> ) 

h  = 

ys 

= 

(J  )h)sin(0)cos(<J>) 

zs 

R 

( li-h)sin(i>) 

L  and  L  in  Rotating  Frame 

The  vectors  L_  and  L_  are  unit  vectors  in  the  earth  centered 
rotating  frame  representing  the  position  and  velocity  of  the  satellite 
relative  to  the  station. 
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V. 


i 

X 

*=^ 

i _ 

)°°b(  0)  _  £32111(0) 

II 

-3l 

h 

- 

(£;>-£ l)ain(0)  t  coa(0) 

u 

R 

eos((J>)co3(e)cos(a)  +  :;in(  ij)nin(e) 

where 


i\  ”  sin((J))eos(e)L03(a) 

(A- 3  ) 

S-2  =  cos($)sin(e) 

(A-4) 

3  =  co3(c)sin(a) 

(A-5) 

( A—  6 ) 


where 

SLX  ::  a[co3(0>3in(^)i5in(d)c:o3(e)  -  cos  ((Os  in  (  O)cos(e)]  + 
e[cos(0)r,.ln(  ;>)sin(e)eos(a)  +  co3(0)eos(<J)co3(e)]  + 
eCsin(  0  )s  i.n(a)  sin(e)  ]  ( A—  7 ) 

l.j  -  dCaiu(0)r,in(ii)  ?in(a)co-'.(e)  +  co3(a)cos(  O)coa(e)]  + 
e[sin(0)Hin(t!>)cos(a)  in(e)  +  r»in(0)cos(  ^)coj(o)]  - 
>[ces( 0)sin(o)sin(a) ]  (A- 8) 

iz  =  -i[co;;(i})«ln(a)sin(e)  1  + 

o[sin($)co:;(o)  -  oos('J>)ccr;(a)3in(e)]  (A-9) 
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Satellite  Position  and  Velocity  in  Rotating  Frame 

Once  the  unit  vectors  L.  and  L.  are  determined,  the  satellite's 
position  and  velocity  in  the  earth-centered  rotating  frame  are 


and 


£r  = 


y 

7. 

- -R 


pL_  +  Rr 


£r  = 


=  pL  t  pL 


(A-10 ) 


( A-ll) 


These  six  components  comprise  the  initial  state  vector  &(0)  in  the 
rotating  system.  The  initial  state  vector  in  the  inertial  system  is 
then  determined  by 


X 

cos(/jjt)  -n.in(wt)  0 

X 

y 

= 

sin(wt)  cos(ut)  0 

y 

7. 

T 

0  0  1 

s 

R 


(A-12) 


and 

/ 
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Appendix  B 

Sta to  Sons itiv.i ty  Matrix  _T 

The  state  sensitivity  matrix  F,  is  used  in  the  computation  of  the 
State  Transition  Matrix  used  in  linearization  schemes.  F  is  defined 
as 


f(x) 

35 


Dx6 

3xj 


3X6 

3x2 


•  •  • 


3*1 

3*1 

Dx, 

3*1 

3x? 

3x6 

3x2 

• 

—  —  •  •  • 

3x2 

3xi 

3x2 

Dx6 

3X6 

3*6 


( B—  1 ) 


and,  when  the  function  f(x)  is  expressed  in  the  inertial  frame,  the 
matrix  F,  neglecting  zero  terms  is 


Fin  "  i'25  -  1’36  =  1 

l'*>  =^[3(v)?  - ’] 

■  7  O 


(13-2) 


( B-3) 


(B-4) 


(B-5) 
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o 

•* 


(n-6) 

(B-7) 

(B-8) 

( B-9 ) 

(B-10) 

(B-ll) 


When  the  function  f(j<)  is  expressed  in  rotating  coordinates ,  the 
F  matrix,  again  neglecting  zero  terms,  is 


F14  "  F25  "  f36  =  1 

( B-12 ) 

F->  =  -rf  [3  (;f  -  *]  *  "2 

( B-13) 

■  ?  (^) 

( B-14 ) 

(B-15) 

f4  5  =  2a) 

( B-16 ) 
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(B-17) 

(B-18) 

(B-19) 

( B-20) 

( B-21) 

( B-22) 


(B-23) 
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Appendix  C 


Measurement  Matrix  M 


The  measurement  matrix  M,  consists  of  the  partial  derivatives  of 
the  measurements  p,  p,  a,  and  e  with  respect  to  the  state  vector  x_.  If 
the  state  vector  is  in  rotating  coordinates,  M  is 


3h( x) " 


(C-l) 


where  M(l),  M(2),  M(3),  and  M(4)  are  row  vectors  and  are 


H(l)  =  ~  [xR-xs,  yR-ys>  sr-zs»  °»  °»  0] 


,,  i  r.  p<xr-xs> 

!)  =  —  xR - , 

P  K  P 


p(yR-ys>  . 

yR  "  n  »  ZR 


p(xR-Xq) 


»  :<r-xs  »  yjj-ys* 


ZR“ZS 


(C-2  ) 


(C-3  ) 


M(3)  =  -  [-XpSin(O)  -  y-j-sinC  )ccs(0),  x-j-cosCG) 

(>’.*  +  y-|; ) 


yTsin(0)sin($),  yTcos($),  0,  0,  0] 


(C-4) 
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M(4)  = 


"1  ~/2 
(xT  +  y^) 


i ( I ) cos ( 0 )  -  ,  sin(0)cos(<}>) 


zr^Ls)  #  sin(l>)  .  SiEpl  ,0,0,0 


(C-5) 


If  the  s cate 


vector  is  in  inertial  coordinates,  the  matrix  M  is 


Hi  =  Mr 


cos(tot)  sin  (ait) 

-sin(wt)  cos  (tot) 


-ojsin(cot)  cos(wt)  0  cos(wt)  sin(iut)  0 

ueos(ut)  -  s in (dit)  0  -sin(ut)  cos(ut)  0 


(C-6) 
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